
#include "kernel/types.hpp"

#include "yocto/math/fit/lsf.hpp"
#include "yocto/sequence/vector.hpp"
#include "yocto/string/vfs-utils.hpp"
#include "yocto/math/io/data-set.hpp"
#include "yocto/ios/icstream.hpp"
#include "yocto/ios/ocstream.hpp"


namespace
{
	
	class Perm
	{
	public:
		const int Z;
		const Real C_in;
		const Real C_out;

		inline Perm( const int __Z, const Real __C_in, const Real __C_out ) :
		Z( __Z ),
		C_in(  __C_in  ),
		C_out( __C_out )
		{
			std::cerr << "Z    =  " << Z << std::endl;
			std::cerr << "C_in =  " << C_in << std::endl;
			std::cerr << "C_out = " << C_out << std::endl;
		}
		
		inline ~Perm() throw() {}
		
		
		Real Compute( const Real zeta, const array<Real> &a )
		{
			assert( a.size() >= 4 );
			const Real zzm   = Z * zeta;
			const Real Psi_m = Psi(zzm);
			const Real SPm   = a[1] * Exp( zeta * a[2] ) + a[3];
			return -Z*Faraday*SPm * Psi_m * ( C_out - C_in * Exp( zzm ) ) + a[4];
		}
		
	private:
		YOCTO_DISABLE_COPY_AND_ASSIGN(Perm);
	};
	
}


int main(int argc, char *argv[])
{
	static const char *progname  = _vfs::get_base_name(argv[0]);
	static const Real  T         = 273 + 25;
	static const Real  F_over_RT = Faraday / ( R*T);
	try
	{
		
		if( argc <= 4 )
			throw exception("%s file.dat Z C_in C_out",progname);
		
		const int  Z     = strtol(argv[2], NULL, 10);
		const Real C_in  = strtod(argv[3], NULL);
		const Real C_out = strtod(argv[4], NULL);
		
		bool use_offset = true;
		bool use_noise  = true;
		for( int i=5; i < argc; ++i )
		{
			const string s = argv[i];
			if( s == "ZOFFSET" )
			{
				std::cerr << "-- " << s << std::endl;
				use_offset = false;
			}
			
			if( s == "ZNOISE" )
			{
				std::cerr << "-- " << s << std::endl;
				use_noise = false;
			}
			
		}
		
		//--- prepare data
		vector<Real> E;
		vector<Real> I;
		vector<Real> Ifit;
		vector<Real> zeta;
		
		//-- load data
		const string source = argv[1];
		const string source_base_name = _vfs::get_base_name( source );
		
		std::cerr << "-- loading data from '" << source << "'" << std::endl;
		{
			data_set<Real> ds;
			ds.use( 1, E); //-- in mV
			ds.use( 2, I); //-- I in pA
			ios::icstream  fp(source);
			ds.load(fp);
		}
		const size_t N = E.size();
		Ifit.make(N,0);
		zeta.make(N,0);
		
		std::cerr << "-- building sample with #data=" << N << std::endl;
		//-- build lsf sample
		for( size_t i=1; i <= N; ++i )
		{
			E[i]   /= 1000.0;           //-- E in V
			zeta[i] = E[i] * F_over_RT; //-- adim
			I[i]   /= 1.0e12;           //-- in A
		}
		
		
		
		fit::sample<Real>      sample( zeta, I, Ifit );
		Perm                   perm(Z,C_in,C_out);
		fit::lsf<Real>::field  field( &perm, & Perm::Compute );
		
		std::cerr << "-- guest shape..." << std::endl;
		{
			const string output = "shape-" + source_base_name;
			std::cerr << "-- saving in '" << output << "'" << std::endl; 
			ios::ocstream fp( output, false );
			for( size_t i=1; i <= N; ++i )
			{
				if( Fabs(E[i]) > 10e-3 )
				{
					const Real zm    = zeta[i];
					const Real zzm   = Z * zm;
					const Real Psi_m = Psi( zzm );
					const Real den   = -Z*Psi_m * ( C_out - C_in * Exp(zzm) );
					fp("%g %g\n", zm, I[i]/den);
				}
			}
		}
		
		
		const size_t nv = 4;
		vector<Real> aorg(nv,0);
		vector<Real> aerr(nv,0);
		vector<bool> used(nv,false);
		
		
		fit::lsf<Real> LSF;
		LSF.ftol = 1e-7;
	
		
		for( size_t i=1; i <= nv; ++i )
		{
			used[i] = true;
			if( !use_offset )
				used[3] = false;
			if( !use_noise )
				used[4] = false;
			
			LSF( sample, field, aorg, used, aerr );
			if( sample.status != fit::success )
			{
				throw exception("Least Square Fit Failure");
			}
			std::cerr << "aorg=" << aorg << std::endl;
			std::cerr << "aerr=" << aerr << std::endl;
		}
		
		std::cerr << std::endl;
		Real &A = aorg[1]; Real &dA = aerr[1];
		Real &B = aorg[2]; //Real &dB = aerr[2];
		Real &C = aorg[3]; Real &dC = aerr[3];
		
		//! A and C is in L/s
		A  *= 1e15; //!< now in microns^3/s
		dA *= 1e15; 
		
		C  *= 1e15;
		dC *= 1e15;
		
		std::cerr << "coef=" << aorg << std::endl;
		std::cerr << "cerr=" << aerr << std::endl;
		{
			const string output = "fit-" + source_base_name;
			std::cerr << "-- saving in '" << output << "'" << std::endl; 
			ios::ocstream fp( output, false );
			fp("#mV pA fit.pA zeta SP\n");
			for( size_t i=1; i <= N;++i )
			{
				E[i]    *= 1000.0;
				I[i]    *= 1e12;
				Ifit[i] *= 1e12;
				fp("%.5e %.5e %.5e %.5e %.5e\n", E[i], I[i], Ifit[i], zeta[i], A* Exp( B * zeta[i] ) + C );
			}
		}
		{
			const string output = "info-" + source_base_name;
			std::cerr << "-- saving in '" << output << "'" << std::endl; 
			ios::ocstream fp( output, false );
			fp( "-- surface*permeability in microns^3/s, from '%s'\n", source_base_name.c_str());
			fp( "\t(%.7f) * exp( (%.7f) * x ) + (%.7f)\n", A, B, C );
		}
		
		return 0;
	}
	catch( const exception &e )
	{
		std::cerr << e.when() << std::endl;
		std::cerr << e.what() << std::endl;
	}
	catch(...)
	{
		std::cerr << "Unhandled Exception!" << std::endl;
	}
	return -1;
	
}